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Abstract 

In recent years the threshold of black hole formation in spherically 
symmetric gravitational collapse has been studied for a variety of mat- 
ter models. In this paper the corresponding issue is investigated for 
a matter model significantly different from those considered so far in 
this context. We study the transition from dispersion to black hole 
formation in the collapse of collisionless matter when the initial data 
is scaled. This is done by means of a numerical code similar to those 
commonly used in plasma physics. The result is that for the initial 
data for which the solutions were computed, most of the matter falls 
into the black hole whenever a black hole is formed. This results in a 
discontinuity in the mass of the black hole at the onset of black hole 
formation. 
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1 Introduction 



The gravitational collapse of a localized concentration of matter to a black 
hole is a central theme in general relativity. Even in the simplest case of 
spherically symmetric collapse much remains to be learned. When no mat- 
ter is present, i.e. in the case of the vacuum Einstein equations, there is 
no collapse in spherical symmetry, because of Birkhoff's theorem. Thus it 
is necessary, in order to obtain information about gravitational collapse by 
studying the spherically symmetric case, to choose a matter model. A simple 
choice, which has provided valuable insights, is that of a massless (real, mini- 
mally coupled) scalar field. A deep mathematical investigation of spherically 
symmetric solutions of the Einstein-scalar field system has been carried out 
by Christodoulou. Some of his results will now be discussed. 

In |J it was shown that for sufficiently small initial data the fields dis- 
perse to infinity at large times. A result for the case of large data was 
proved in |J. For any solution it is possible to define a number M with the 
property that in the region r > 2M (where r is the area radius) the solution 
approaches the Schwarzschild solution of mass M at large times. The phys- 
ical interpretation is that the system has collapsed to form a black hole of 
mass M. Of course M = in the case that the field disperses. The dichotomy 
between dispersion and black hole formation still leaves the question unan- 
swered, which data result in which of these two outcomes, except for small 
data. In |7| Christodoulou gave a sufficient condition on initial data to ensure 
that M > for the corresponding solution, although this criterion is not very 
practical. 

For large data the above results only describe the structure of the solution 
at sufficiently large radius and the internal structure is left open. In || it was 
proved that there exist data leading to the formation of a naked singularity. 
The negative implications of this for the cosmic censorship hypothesis are 
limited by the fact || that this behaviour is unstable in the class of spherically 
symmetric initial data. 

An important new element was brought into the study of the Einstein- 
scalar field and gravitational collapse in general by the numerical work of 
Choptuik ||. He took a fixed initial datum for the scalar field (which due 
to spherical symmetry determines initial values for the gravitational field) 
and scaled it by an arbitrary constant factor. This gives rise to a family of 
initial data depending on a real parameter A. For all initial data used in the 
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computations the same picture emerged. For values of A corresponding to 
small data the field dispersed, in agreement with the rigorous results. For 
large initial data a black hole was formed. Under these circumstances we 
can define a critical parameter A* as the lower limit of those values of A 
for which a black hole is formed. If the mass M of the black hole formed is 
plotted as a function of A it is found that M(A) is continuous. In particular 
lim^^ t M(A) = so that black holes of arbitrarily small mass can be formed 
within the given one-parameter family. Up to now it has not been possible 
to confirm this behaviour by rigorous mathematical arguments. Choptuik's 
results contain much more detailed statements than the one just mentioned, 
but it is the one which will be important in the sequel. 

The nature of the boundary between dispersion and black hole formation 
in gravitational collapse is now a very active research area. Most (but not 
all) of the work has been concerned with spherical symmetry and varying 
matter models and has relied essentially on numerical computations. For 
a recent review of the field see fTTJ. In a study of critical collapse for the 



Einstein- Yang-Mills equations, Choptuik, Chmaj and Bizoh |4j] have found 
both cases where \im.A-^A»M(A) — and ]im.A->A*,A>A* M(A) > 0. They call 
these cases type II and type I respectively, the terminology being motivated 
by an analogy to phase transitions in statistical mechanics. They relate 
type I behaviour in this system to the existence of the Bartnik-McKinnon 
solutions H], |23[], which are static. (The models which had been considered 
previously, and which showed exclusively type II behaviour, admit no regular 
static solutions.) 

Almost all the matter models which have been considered up to now 
in the context of critical collapse are field theoretic in nature rather than 
phenomenological. The one exception is a perfect fluid with linear equation 
of state, which shows type II behaviour. A phenomenological matter model 
whose collapse is of interest is collisionless matter described by the Vlasov 
equation. As will be described next, the spherical collapse of collisionless 
matter has been studied both analytically and numerically but the known 
results for this type of matter say little about the nature of critical collapse. 
The purpose of this paper is to begin the numerical investigation of critical 
collapse of collisionless matter. 

In fTM it was shown that for sufficiently small initial data for collisionless 
matter the matter disperses to infinity at large times. Thus the analogue of 
Christodoulou's small data result holds for collisionless matter. The solutions 
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are geodesically complete. Unfortunately, no analogue of his large data result 
is known is this case. The only result in that direction which has been proved 
is that there do exist initial data which develop singularities W71 . The proof 



proceeds by demonstrating the existence of initial data which contain trapped 
surfaces and applying the Penrose singularity theorem. A different kind 
of large data result, which is relevant to the numerical calculations of this 
paper, is that if data given on a hypersurface of constant Schwarzschild time 
gives rise to a solution which develops a singularity after a finite amount of 
Schwarzschild time, then the first singularity occurs at the centre of symmetry 



An analogous result where Schwarzschild time is replaced by maximal 



slicing has also been proved [|T| 



It seems plausible that in those solutions which develop singularities black 
holes are usually formed. However there are no mathematical results on 
this and the convincing evidence for the formation of black holes is purely 
numerical. This is part of a large body of work due to Shapiro, Teukolsky 
and collaborators which goes far beyond the spherically symmetric case. Here 
only those results will be discussed which are directly relevant to this paper. 



In particular only results for spherical symmetry are covered. In f20| , |21| , yH 
the collapse of various configurations of collisionless matter to a black hole has 
been computed numerically. The relative merits for this problem of different 
choices of time coordinate (polar or maximal) and radial coordinate (area or 
isotropic) are discussed in [[12], p~3|j . 

In this paper we carry out an experiment analogous to that of Chop- 
tuik for spherically symmetric collapse of collisionless matter. We use 
Schwarzschild coordinates (i.e. polar slicing and area radius). Starting with 
a suitable smooth function fo of compact support as initial datum for the 
distribution function we consider the scaled data Afo, where A is a positive 
constant and compute the corresponding time evolution numerically. For all 
data for which we tried the experiment, the results were qualitatively sim- 
ilar. When A is sufficiently small the matter disperses, in agreement with 
the analytic theory. This happens up to some value A* of A. For A > A* 
the observed behaviour indicates the formation of a black hole. The lapse 
function develops an abrupt step at a certain radius r(A). This step remains 
at the same radius but gets deeper and deeper. We interpret this as the 
signature of a black hole with mass M(A) =r(A)/2. If M(A) is plotted as 
a function of A it is found that the limiting value of M(A) as A approaches 
A* from above is strictly positive. Thus we find behaviour of type I in the 



4 



terminology of ||]. We never find any signs of singularity formation for any 
value of A and this is consistent with the standard picture where the only 
singularities formed are those of black hole type and they are avoided by a 
Schwarzschild time coordinate. 

As a check on the interpretation of the numerical solutions as describing 
collapse to a black hole, radial null geodesies were computed. The results 
agreed well with the expected picture. Radial null geodesies starting at the 
centre at early times escape to large values of r. Those starting after a 
certain time T\ remain within a finite radius. The limit of this radius as t 
tends to T\ from above is equal to r(A) as computed above. Thus we obtain 
a consistent picture with a black hole whose event horizon is generated by 
the null geodesies starting at the centre at time T±. 

To have a better picture of what is happening in the collapse of near 
critical initial data it is useful to consider how the mass M(A) of the black 
hole formed depends on the ADM mass of the initial configuration. For 
parameter values A well above A* these quantities are almost equal. In other 
words, essentially all the matter falls into the black hole. As the critical 
parameter value is approached from above some of the matter does not fall 
into the black hole and in some of the cases which were computed escapes 
to infinity. However, in all the cases which were computed the mass of the 
black hole is more than 90% of the total ADM mass of the configuration. 
Thus whenever a black hole is formed almost all the matter falls into it and 
the mass gap is a reflection of this. The picture in the only other case of a 
phenomenological matter model in which critical collapse has been studied up 
to now, namely a perfect fluid with linear equation of state, is very different. 
In the case of a radiation fluid (equation of state p=(l/3)p) the slightly 
supercritical collapse can be described as follows flnj. The matter divides 



almost completely into two parts, separated by a near vacuum region. The 
outer part of the matter, which contains almost all the mass, escapes to 
infinity. The inner part, which contains only a very small amount of mass, 
collapses to form a (small) black hole. As the critical parameter is approached 
this situation becomes more and more extreme and the black hole mass tends 
to zero. 

The numerical code used is based on a numerical scheme for the corre- 
sponding Newtonian problem described in f!9]| . It incorporates less refined 
features than the codes of Shapiro and Teukolsky but seems to be quite suf- 
ficient for the present task. It does have the advantage that an analogous 



5 



Newtonian code has been proved to be convergent [19| and it seems reason- 
able to hope that this proof could be extended so as to obtain results on the 
convergence of the method used in this paper in the future. 

These results add collisionless matter to the class of matter models for 
which something is known about critical collapse. Clearly it is desirable to 
examine other types of initial data so as to discover the prevalence of the 
type of behaviour found here or of others which have not yet been seen. It 
is our hope that these numerical investigations can also help to further the 
mathematical study of collisionless matter in general relativity by providing 
pictures of what is happening which can suggest which theorems one should 
try to prove and by what means. Mathematical investigations of partial 
differential equations often proceed by estimating the growth rates of certain 
quantities and it is useful to have an idea of the expected growth rates on the 
basis of numerical computations. The fact that we observe no singularities 
in the numerical computations can be seen as evidence that the weak cosmic 
censorship conjecture is true for collisionless matter. Indeed it may even be 
true in a stronger version than in the case of the massless scalar field. There 
may be no naked singularities formed for any regular initial data rather than 
just for generic initial data. This speculation is based on the fact that the 
naked singularities which occur in scalar field collapse appear to be associated 
with the existence of type II critical collapse. 

The paper proceeds as follows: In the next section we formulate the 
Vlasov-Einstein system, first in general coordinates, and then in coordinates 
adapted to the spherically symmetric, asymptotically flat situation that we 
want to study. In Section 3 we describe the code we are using and test it on a 
steady state. In Section 4 we present the results of the numerical simulations. 

2 Formulation of the spherically symmetric 
Vlasov-Einstein system 

In the present paper the matter model is a collisionless gas as described by 
the Vlasov or Liouville equation. Coupling this equation self-consistently to 
the Einstein field equations results in the Vlasov-Einstein system, which we 
first write down in general coordinates on the tangent bundle TM of the 
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spacetime manifold M: 

p a d xa f-r«y P w pa f = o, 

G af3 = 8nT al3 , 
J m 

Here / is the number density of the particles on phase-space, and G a/3 de- 
note the Christoffel symbols and the Einstein tensor obtained from the space- 
time metric g a p, \g\ denotes its determinant, T alS is the energy-momentum 
tensor generated by /, x a are coordinates on M, (x a ,pP) the corresponding 
coordinates on the tangent bundle TM, Greek indices run from to 3, and 

rn=\g afi p a pP\ 1/2 

is the rest mass of a particle at the corresponding phase-space point. We 
assume that all particles have rest mass 1 and move forward in time so that 
the distribution function / lives on the mass shell 

PM = {g aP p a T P = -l, p°>0}. 

We consider this system in the asymptotically flat and spherically symmetric 
case and use Schwarzschild coordinates to coordinatize the spacetime mani- 
fold. The metric takes the form 

where t e H, r > 0, 6 e [0, 7r] , G [0, 27r] . Asymptotic flatness is then expressed 
as the boundary condition 

lim \(t,r) = lim u(t,r) = 0. 

r^oo r— too 

We also require a regular centre, which is guaranteed by the boundary con- 
dition 

A(t,0) = 0. 

To write the Vlasov equation we use the corresponding Cartesian coordinates 

x— (rsin0cos^,rsin0sin^,rcos0) 



7 



as spatial and 



v * = p* + (e x -i)^-—, a = l,2,3 



as momentum coordinates. The Vlasov-Einstein system then takes the form 
3 m-a v _ a f_( \^ , m-a, ,/ /- i .<A x 



d t f + e»~ A L _ 2 -^/-(A— + e^-VV 1 + H a )r-^/ = 0, (2-1) 



e~ 2A (2rA'-l) + l = 8vrr 2 p, (2.2) 
e- 2A (2?y + 1) - 1 = 87rr 2 p, (2.3) 



p(t,r) = p(t,x) = y yi + H2/(t s x,v)di;, (2.4) 
p t,r = p t,x = / / t,x,^ . 2.5 

7 v r ^ \A + H 2 

Here x,wGlR 3 , r=|x|, • denotes the Euclidean dot product in 1R 3 , and A = 
dX/dt and // = d\ijdr. f is assumed to be spherically symmetric in the sense 
that 

/(<,&,!;) = /(*, Ac, AgSO(3). 

Eqns. (|2.2f) and ( |2.3|) are the 00- and 11-components of the field equations; 
it can be shown that for a solution of the above system also the remaining 
nontrivial components of the field equations hold. We state the 01-component 
explicitly, since it is used in our numerical scheme: 

A = -4vrre A+/i j, (2.6) 

where 

j(t,r) =j(t,x) = [ —f(t,x,v)dv. (2.7) 



This form of the spherically symmetric Vlasov-Einstein system is convenient 
for analytical work and has been used in p4| , |16|j . For numerical work one 



wishes to use the symmetry explicitly in the Vlasov equation in order to 
reduce the latter's dimension. One set of independent variables to use is 

-i x ' v 

r=\x\, u—\v\, a = cos . 

ru 
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However, an equivalent and more convenient set of variables is 

r=|x|, w = — — , L— \x\ 2 \v | 2 — (x-v ) 2 = \x x t>| 2 , (2.8) 
r 

particularly because L is constant along the characteristics of ( |2.1| ); note that 

2 2 , L 
u =w +—. 

In these variables the Vlasov equation for f = f(t,r,w,L) becomes 

d t f + e »- x ^JL=d r f - (\w + ^Vv/I+^-e^ fl^ / = 0. 

Vl+M \ rVl+M / 

(2.9) 

The field equations remain unaffected, and the source terms (|2.4j ), ( |2.5|) , (|2.7|) 
can be rewritten in terms of (r,«,a) or (r,w,L). 

We now mention some results that have been established for the spheri- 
cally symmetric Vlasov-Einstein system. Throughout we consider as initial 

o 

condition a spherically symmetric, nonnegative function / which as a function 
of x and v is continuously differentiable, has compact support and satisfies 
the inequality 

/ / Vl + v 2 f{x,v)dvdx<-, r>0, (2.10) 

J\x\<rJ 2 

which means that the initial hypersurface does not contain a trapped surface. 
In |14| it was shown that for such an initial condition there exists a unique, 

o 

continuously differentiable solution / with /(0) = /, which exists with respect 
to Schwarzschild time on some right maximal interval [0,T[. If the solution 
blows up in finite time, i. e., if T < oo then p(t) becomes unbounded as t — > T~ . 
Actually, as shown in p in this case has to become unbounded at the 
centre r = 0, i. e., if any singularity evolves, the first one must be at the centre. 
This rules out singularities of shell crossing type, which can be a nuisance 
in other matter models, e. g. dust. If the initial datum is uniformly small 
the resulting solution is global in the sense that the spacetime is geodesically 
complete and the components of the energy momentum tensor as well as 
metric quantities decay with certain algebraic rates in t. 
Let us denote 

m(t,r) =4vr / s 2 p(t,s)ds. (2.11) 
Jo 
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Then m(t, oo) is a conserved quantity, the ADM mass of the system. Using 
m(t,r) the field equations (|2~2| ) and (|2.3|) yield 



-2A(t,r) _ 



2m(t,r) 



//(t,r) 



3 2A(t,r) 



m(t,r) 



■4nrp(t,r) 



(2.12) 
(2.13) 



also 



A'(t,r) = ^;(e 2A (87rr 2 p(t,r)-l)-l) 



„2A 



^M+4vrrp(t,r)V (2.14) 



note that the right hand side of ( 2.12 ) is positive at t = by the assumption 
( |2.10| ). A further quantity which is conserved by the system is the total 
number of particles 

e x f(t,x 1 v)dvdx. (2.15) 



3 Description and Testing of the Code 

o 

Let us consider an initial condition, f(x,v), which is spherically symmetric, 
satisfies the condition fl2.1U| ), and vanishes outside the set (r,u,a) G [R ,Ri] x 
[C/ ,^i] x [a ,ai]. We will approximate the solution using a particle method. 
For a thorough treatment of particle methods in the context of plasma physics 
see [Qj. To generate the particles we take integers N r ,N u ,N a and define 

Ri-Ro Ui-Uo «i-« 

r i = R +(i- ^ Ar, Uj = U + (j - Au, a k = a +(k- ^ Aa, 

fijk = f( r ii u j, a k) 47rrfAr 27ra| Au sina fc Aa, 
r °ijk = r h Wi jk = UjCOsa k , L ijk = (r i u j sma k ) 2 . 

From these, approximations are made of the quantities p, p, j, and m defined 
in (|2.4|), ( |2.5|) , ( |2.7| ), and ( |2.11| ) at the grid points nAr. The equation 

1 / Tfl \ 
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which is obtained from ( |2.12| ) and ( |2.13| ), together with ^-+0 as r->oo is 
used to compute \i on this grid. Note that for r outside the support, p = and 
m = constant and this equation is explicitly integrable. A is computed using 
(|2.12|) , and similarly A and A' are computed using ( p^6|) and ( |2.14|) . Letting 
D denote differentiation along a characteristic of the Vlasov equation (|2.9| ) 
we have 

Dr = e»- x " 



Dw = -wX - e^Vl + uV + ■ 



recall that 



r 3 Vl + u 2 
DL = 0: 



2 2 , L 



We interpolate p, j, p, p, A etc. to particle locations and use these equations 
to define r\- k and w}^ k using a simple Euler time stepping method. Here r\- k 
denotes an approximation of the radius of the characteristic at time At with 
r = r i) u = Uj, a = ctk at time 0. 
We also have the equation 



( At) - 1 (fl jk - = -% k (A + w\'e^ x (l + 



which represents the time evolution of a volume element along a character- 
istic. One time step is now complete. 

To test the code a steady state solution was generated. Following |L5] we 
take 

f(x,v) = <f>(E) 

where 

E = e^^Jl + \v\ 2 
is the particle energy. For simplicity 



1 , E < E 
, E > E 



4(E): 

with Eq > was taken. Then from |15| , 

p{r)=9^{K r )) 
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and 

p{r) = hf(ji(r)) 
where _ 

0(ee")eVe 2 -lrfe 

and 

^H = y/°°0(6e«)( e 2 -l) 3/2 rf6. 
Substituting into (R7TJ), /x must satisfy 

1 /47T /' r \ 

" V) ° l-HW))^ ^l ^W.))*+4«r»*Wr))) (3.2) 

and /j— >0 as r^oo. This was solved using a shooting method. Note that 
for /i>lni? ( P. 2|) reduces to 

,, . r~ 2 m(oo) 

/x (r) = : — - — -, 

1 — 2r i m(oo) 

which is explicitly integrable. Thus, it was only necessary to solve (|3.2|) on a 
bounded domain. Eq = 0.9 proved to be a convenient choice. The resulting 
steady state has mass 3.36 10~ 2 and support contained in 0<r<0.36. The 
radial component of v ranges from —0.64 to 0.64. The maximal values of \i 
and A are 0.296 and 0.132 with p a decreasing function of r. 

It was found that taking iV r = 40, iV u = 10, and N a = 10 produced good 
results. This resulted in 2550 particles (less than 40 x 10 x 10 since the sup- 
port of / is not rectangular). At time zero the maximal errors in m and \i 
(maximum over r) were 1.63 10 -4 and 2.84 10 -3 respectively. Dividing by 
the maximal values of m and \i (that is, by 3.36 10 -2 and 2.96 10 _1 ) we find 
the maximal errors in m and \x are 0.49% and 0.96% at t — 0. 

The errors in m and /i were computed at time t = 10 using different time 
steps. Percentage errors were computed as described above with the results: 
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At 


error in m 


error in fi 


1 

4000 


6.2% 


5.9% 


1 

8000 


3.4% 


2.9% 


1 

16000 


2.1% 


1.3% 



Thus the particle code tracks the steady state reasonably well, although a 
rather small time step seems to be needed. We attribute this, at least in part, 
to numerical difficulties in tracking the motion of particles near r = 0. Note 
that for this steady state the density is largest at r = 0. For other solutions 
with zero density near r = 0, the time step was taken larger without significant 
change in the results. 



4 Results of Simulation 

In this section we consider initial data 

o 

f(x,v)=Af (x,v) 
with fo fixed and vary A. As a first example we take 

f (x,v)= [50,000(2.2 -r)(r -2)(10.2-w)(w-10)(3.1-a)(a-2.9)] 2 

for 2<r<2.2, 10<w<10.2, 2.9<a<3.1 and f o (x,v) = otherwise. Thus 
the mass is initially concentrated between r = 2 and 2.2 and is moving inward 
rapidly. In most of the following simulations the support of / is divided into 
40 by 20 by 20 cells (40 in r) resulting in 16000 particles and the time step 
is 0.005. The cases where this is not so will be pointed out. 
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In figures 1 and 2, A is 0.69. In figure 1 the enclosed mass m(t,r) defined 
in ( j2.11| ) is plotted at times t = 0, 2, 4, and 8. 



Figure 1: Enclosed Mass, Subcritical 
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Figure 2 shows fi at the same times. The solid curve is t = 8 in both. At 
time 8 every particle is moving outward with momentum greater than 3.37 
and has position r > 1.48. In figure 1 we see that the mass has fallen inward 
and reversed direction, almost returning to its starting position at time 8. 
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Comparing figures 1 and 2 we see that at t = 4 the mass is near r = 
and that attains relatively large values. By t = 8 fi has dropped to near 
its starting values. As t grows the particles continue outward and disperse, 
consistent with the small data result [|14| . 



Figure 2: Mu, Subcritical 
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In figures 3 and 4, A is 0.75. Again m and fi are plotted at times t = 0, 2, 4, 
and 8 with the solid curve representing t = 8. For times 0, 2, and 4 figure 3 
resembles figure 1 , but for t = 8 we see in figure 3 that the mass has not moved 
back out. Rather it is centered near 0.21, and from figure 4 fi has formed 
an abrupt transition near 0.235. At time 8 the maximal outward momentum 
is about 9 while the maximal inward momentum is 189 and the positions 
satisfy 0.13 < r < 0.24. Examination of j{t,r) defined in (|2.7| ) reveals that the 
mass flux is almost entirely inward. 
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Figure 3: Enclosed Mass, Supercritical 
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The abrupt transition in /i occurs at r = 0.235 and A has a maximum at this 
same r. A forms a cusp at its maximum. At time 16 these features have not 
moved although the maximal value of has grown. 
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Figure 4: Mu, Supercritical 
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Recall that the quantities 



and 



J J e x fdvdx 

J J \Jl + \v\ 2 fdvdx 

are both conserved by the exact time evolution. In the above runs neither 
quantity varied by more than 1.7% of its initial value. A run was made with 
80x40x40 particles, At = 0.0025, A = 0.75, and final time 8. The resulting 
graphs of m and /i with the finer resolution are qualitatively very similar to 
figures 3 and 4, except that the transitions are slightly more abrupt and the 
maximal value of \fj,\ is increased by about 7%. 

When A was taken larger than 0.75 the results are similar. For 0.70 < 
A< 0.74 a similar structure formed with a stationary abrupt transition in //, 
but a small amount of mass escaped. Thus it seems that the critical value of 
A for this choice of fo is A* ^ 0.70. 

For A > 0.70 the radius of the maximum of A and the radius of the abrupt 
transition in \i are nearly the same. Thus we have computed the radius r(A) 
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where the maximal value (over r and t) of 



A: 



1, / 2m\ 
— n 1 
2 V 



r / 



occurred. So this is the value of r where the condition ( [2.10| ) is most nearly 
violated. For A > 0.70 the maximum A was attained at the largest time, for 
A < 0.69 it was attained earlier. The results are graphed in figure 5. 



Figure 5: Radius vs A 
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Since the ADM mass of the configuration depends linearly on A and for 
A > 0.70 nearly all the mass is captured in the black hole, the mass of the 
black hole, M(A) = r(A)/2, depends nearly linearly on A for A > 0.70. The 
values for A < 0.70 are plotted as a dotted curve. We note that 

lim M(A) = - lim r(A)w0.11, 

A~,A+ 2 A^A+ 

which indicates type I behaviour as explained in the introduction. Since 
Ar = 0.005, the discontinuity in M(A) at A* is significant. 
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Next we consider 



/o(^) = 0.l(l-r 2 ) 2 (l- M 2 ) : 



for r<l and u<l and fo(x,v) = otherwise. Similar to figure 5, figure 6 
shows the radius where the maximal value of A occurs as a function of A. 
For A> 1.6 the maximum was attained at the largest time. For this initial 
condition a smaller time step of 0.00125 was used. 

Figure 6: Radius vs A 
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Similarly we consider 

/o(^,^) = 0.1(3-r) 2 (2-r) 2 (l-r) 2 (l- M 2 ) 2 

for 1 < r < 3 and u < 1 and f (x,v) = otherwise. Figure 7 plots radius versus 
A as figures 5 and 6 did. For A > 0.76 the maximum of A occurred at the 
largest time. For A < 0.75 the maximum occurred at times near zero, hence 
the nearly constant values of r for A< 0.75. 
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Figure 7: Radius vs A 
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In figures 5, 6 and 7 the final time was taken large enough that increasing it 
produced only minor changes. In each case we see that the radius at which A 
is largest and the step in the lapse function e 2 ^ forms remains bounded away 
from zero. 
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